Differential expression analysis

The following plots show the 20 most differentially expressed genes for each cluster. Note that the comparisons are within a certain cell-types only (i.e. the CD4+ sub-clusters are compared against other CD4+ subclusters, but not against CD8+ and NK cells).

GO-term enrichment analysis

Differntial expression analysis

We compare gene expression of HPV+ vs HPV-. We perform the comparision * for all cell-types together (=“overall”) * for each major cell type individually

Note that for CD8+ T cells, the major differences (TRBV23, TRBV27) are mainly driven by the cluster of CD8+ T cells that are only abundant in a single patient (CD8_2).

Red-ish genes are over-represented in HPV16+, blue ones under-represented.

Differential expression by bulk

The enrichment of HSPA1A and CCL4 holds even in the bulk model. Yet, it is not statistically significant.

NKRs in single cell data

nkrs = [
    "TNFSF10",  #TRAIL
    "FCGR3A",   #CD16
    "NCR3", #NKp30a,b
    "KLRC2", #NKG2C
    "KLRK1", #NKG2D
    "CD244", #2B4
    "CD226", #DNAM-1
    "KLRC1", #NKG2A
    "KIR3DL1", #KIR
    "TNFRSF9", #CD137/41BB
    "TNFRSF4", #OX40
    "CD27", #CD27
]

adata_cd8 = adata[adata.obs["cell_type"] == "T cell CD8+", :]
sc.pp.highly_variable_genes(adata_cd8, flavor="cell_ranger", n_top_genes=2000)

## extracting highly variable genes
##     finished
## added
##     'highly_variable', boolean vector (adata.var)
##     'means', float vector (adata.var)
##     'dispersions', float vector (adata.var)
##     'dispersions_norm', float vector (adata.var)
## Trying to set attribute `.var` of view, making a copy.

sc.pp.neighbors(adata_cd8, n_neighbors=10)

## computing neighbors
##     using 'X_pca' with n_pcs = 50
##     finished

sc.tl.umap(adata_cd8)

## computing UMAP
##     finished

sc.pl.umap(adata_cd8, color=["hpv_status", "ir_status"] + nkrs)

## /data/scratch/sturm/scratch/projects/2020/kortekaas2020_paper/work/conda/run_notebook-b99176dfb05e0522988294512373169f/lib/R/library/reticulate/python/rpytools/call.py:13: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
##   res = rpycall.call_r_function(f, *args, **kwargs)

sc.pl.dotplot(adata_cd8, nkrs, groupby="hpv_status")

## GridSpec(2, 5, height_ratios=[0, 10.5], width_ratios=[4.199999999999999, 0, 0.2, 0.5, 0.25])

sc.pl.dotplot(adata_cd8, nkrs, groupby="ir_status")

## GridSpec(2, 5, height_ratios=[0, 10.5], width_ratios=[4.199999999999999, 0, 0.2, 0.5, 0.25])

Comparing all 13 tumors: what T cell subgroups can we find?

We independently clustered CD8+ T, CD4+ T and NK cells using unsupervised clustering with the Leiden-Algorithm. In total, we identified 10 clusters, of which there are 3 NK cell subtypes, 3 CD4+ T cell subtypes and 4 CD8+ T cell subtypes. It is an inherent limitation of unsupervised clustering that the “true” number of clusters is not known a priori and the number of detected clusters is somewhat arbitrary. If a more detailed view of subtypes/states should be required, the resolution and thereby the number of clusters can be increased.

The following figure shows a UMAP plot of all T and NK cells colored by cluster.

We first asked, if the clusters are equally distributed across patients. The following plots show

• For each cluster, the fraction of cells originating from a certain patient
• For each patient, the fraction of cells falling into a certain cluster

We observe that some clusters, in particular CD8_2 and NK_1 are mostly driven by a single patient. While we cannot rule out entirely that this is driven by batch effects, it seems unlikely because:

• These patients are outliers only for the given cell-types (e.g. CD8+) but integrate well for the other cell types (CD4+ and NK).
• While the cluster seems unique for the patient, the respective patient-sample also contains T/NK cells that fall into other clusters.

Can we see differences between HPV- and HPV+ tumors

The following plot shows the fraction of cells attributed to a certain cluster for HPV+ and HPV- patients respectively. The black dot refers to the median, the black vertical line expands to the median absolute deviation (MAD). The p-value has been calculated using a Wilcoxon-Mann-Whitney test.

Note that this analysis is not ideal, as it depends on the number of cells detected per patient. This analysis assumes that the number of detected NK/T cells per patient is proportional to the true NK/T cell inflitration. This is reasonable, as we obtained the cells from a dataset containing also other cell types – we can assume that the NK/T cell infiltration is higher if they make up a larger fraction of the total cells.

On the other hand, the number of detected cells per sample also depends on the efficacy of FACS pre-sorting and the library preparation. However, these should be random processes and not related to the patient group.

We performed differntial gene expression analysis between HPV+ and HPV- patients. For each patient, we aggregated all cells of a certain type to an “artificial bulk” RNA-seq sample and compared between those samples.

(I do not perform this DE-analysis on the single cell level, as the HPV-status is confounded by the patient. If I performed this analysis on the single-cell level, we would not get reliable p-values as the DE-analysis tool treats each cell as a sample. Even if a change was driven by a single patient only, the expression difference would appear as highly statistically significant)

None of the genes is statistically significantly differentially expressed after multiple-testing correction. The most differentially expressed genes are shown in the following plot. Due to the small sample size (n=13 patients) they might still be interesting candidates for follow-up.

Can we see differences between HPV16 immune response-negative and - positive tumors in HPV16+ tumors?

The following plot shows the fraction of cells attributed to a certain cluster for IR+ vs IR- patients. The black dot refers to the median, the black vertical line expands to the median absolute deviation (MAD). The p-value has been calculated using a Wilcoxon-Mann-Whitney test.

Again, note that this analysis is not ideal due to the reasons stated above.

We performed differntial gene expression analysis between IR+ and IR- patients. For each patient, we aggregated all cells of a certain type to an “artificial bulk” RNA-seq sample and compared between those samples.

None of the genes is statistically significantly differentially expressed after multiple-testing correction. The most differentially expressed genes are shown in the following plot. Due to the small sample size (10 patients) they might still be interesting candidates for follow-up.

Can we find NK receptors on CD8 T cells? Differences in tumor types?

I retrieved a list of NK receptors from https://www.frontiersin.org/articles/10.3389/fimmu.2015.00601/full (see figure below) and mapped them to gene symbols manually using genecards.org.

NKR figure from the review paper:

GO term enrichment: Biological process

The plots show the 20 most enriched GO-terms from the biological-process ontology. The further right a point, the higher the score of the corresponding term.

GO term enrichment: Molecular function

The plots show the 20 most enriched GO-terms from the molecular function ontology. The further right a point, the higher the score of the corresponding term. If fewer terms are shown, there were fewer terms that met the statistical significance threshold.